Podsumowanie

Średnia długość życia wydaje się prostą do przewidzenia zmienną dla krajów, posiadających dane statystyczne. Klasyfikacja na poziomie 0.95 R2 jest dostatecznie wysoka, aby wykorzystywać dane statystyczne do zwiększania średniej życia w krajach – np. poprzez inwestowanie pięniędzy bądź czasu w czynniki, które najbardziej na to wpływają – są to między innymi poziom nauczania i ilość zgonów spowodowanych wirusem HIV.

1. Wyliczanie bibliotek

Użyte bibliotek przedstawiam poniżej. Są to biblioteki do operacji na macierzach, statystyczne, do uczenia maszynowego oraz do przedstawiania wykresów.

usePackage <- function(p) {
    if (!is.element(p, installed.packages()[,1]))
        install.packages(p, dep = TRUE, verbose = FALSE)
    require(p, character.only = TRUE, quietly = TRUE)
}

packageList <- c("dplyr", "matrixStats", "tidyr", "caret", "ggcorrplot", "plotly")

print(packageList)
## [1] "dplyr"       "matrixStats" "tidyr"       "caret"       "ggcorrplot" 
## [6] "plotly"

Instalacja bibliotek obdywa się za pomocą skryptu, z wyłączeniem logu.

sapply(packageList, usePackage)
##       dplyr matrixStats       tidyr       caret  ggcorrplot      plotly 
##        TRUE        TRUE        TRUE        TRUE        TRUE        TRUE

2. Powtarzalność obliczeń

Powtarzalność obliczeń jest zapewniona najprostszym sposobem – wywoływaniem funkcji, która ustawia ten sam seed dla generatora liczb losowych.

setConsistency <- function() {
  set.seed(42)
}

3. Ładowanie danych

Ładowanie danych polega na pobraniu ich z internetu, a następnie załadowanie jako tabeli.

loadData <- function(urlPath, localFile = 'life_expectancy.tsv') {
  remote_data_file <- download.file(urlPath, destfile = localFile, quiet= TRUE)
  return(read.csv(TMP_LOCAL_FILE, quote='"'))
}
FILE_URL <- 'http://www.cs.put.poznan.pl/alabijak/emd/projekt/Life_Expectancy_Data.csv'
TMP_LOCAL_FILE <- 'life_expectancy.csv'

data <- loadData(FILE_URL, TMP_LOCAL_FILE)

4. Przetwarzanie brakujących danych

Podczas wstępnego przetwarzania danych zmieniane są nazwy kolumn, które są niejednolite. Dodatkowo, spożycie alkoholu oraz procent rządowych wydatków na ochronę zostały wyzerowane, jeśli były “NA”, a rok był 2015. Te przykłady nie były usuwane, gdyż każdy przykład nie posiadał tych danych, co może sugerować po prostu na brak danych z tych lat, a większosza ilość danych pomaga klasyfikatorom.

Następnie,przykłady posiadające “NA” w jakimkolwiek atrybucie zostały usunięte, a także utworzony został dodatkowy zbiór posiadający jedynie atrybuty liczbowe.

fixData <- function(df) {
  # Rename data
  fixedData <- rename(df,
     Life.Expectancy = Life.expectancy,
     Infant.Deaths = infant.deaths, 
     Percentage.Expenditure = percentage.expenditure,
     Under.Five.Deaths = under.five.deaths,
     Total.Expenditure = Total.expenditure,
     Thinness.Years.10.19 = thinness..1.19.years,
     Thinness.Years.5.9 = thinness.5.9.years,
     Income.Resources.Composition = Income.composition.of.resources)
  
  fixedData <- within(fixedData, Alcohol[is.na(Alcohol) & Year == 2015] <- 0)
  fixedData <- within(fixedData, Total.Expenditure[is.na(Total.Expenditure) & Year == 2015] <- 0)
  
  return(fixedData)
}

cleanData <- function (df) {
  return(na.omit(df))
}
  
removeChrColumns <- function(df) {
  library(dplyr)
  return(df %>% 
    select_if(~!is.character(.)))
}

fixedData <- fixData(data)
cleanedData <- cleanData(fixedData)
numericData <- cleanData(removeChrColumns(cleanedData))

5. Rozmiar zbioru i statystyki

Rozmiar danych przedstawiony jest poniżej.

cat(paste('Zbiór danych zawiera', nrow(cleanedData), 'przykładów.\n'))
## Zbiór danych zawiera 1777 przykladów.
cat(paste('Każdy z przykładow posiada', ncol(cleanedData), 'atrybutów.\n'))
## Kazdy z przykladow posiada 22 atrybutów.

Do podstawowych statystyk zostało wybrane tylko kilka atrybutów. Jak możemy zaobserwować, średnia długość życia wacha się od 44 do 89 lat, ze średnią na poziomie 69 lat oraz medianą 71 lat. Jest to pozytywna wiadomość, gdyż oznacza, iż krótka długość życia występuje w niewielu krajach. Widać także bardzo duże dysproporcje w liczbie zgonów dzieci oraz na wirusa HIV.

library(dplyr)
library(tidyr)

cleanedData %>%
  tibble::as_tibble() %>% 
  select(Life.Expectancy, Schooling, Infant.Deaths, HIV.AIDS) %>%
  summarise(across(everything(), list(mean = mean, median = median, min = min, max = max), .names = "{.col}_{.fn}")) %>%
  gather(variable, value) %>%
  separate(variable, c("var", "funkcja"), sep = "\\_") %>%
  spread(var, value) %>%
  relocate(Life.Expectancy, .after = funkcja)
## # A tibble: 4 x 5
##   funkcja Life.Expectancy HIV.AIDS Infant.Deaths Schooling
##   <chr>             <dbl>    <dbl>         <dbl>     <dbl>
## 1 max                89      50.6         1600        20.7
## 2 mean               69.4     1.90          32.2      12.2
## 3 median             71.7     0.1            3        12.3
## 4 min                44       0.1            0         4.2

6. Szczegółowa analiza atrybutów

Przedstawione są histogramy wszystkich atrybutów. Możemy na nich zaobserwować niepokojące ekstrema – wygląda to tak, jakby niektóre z krajów żyły w ciągłej biedzie, a niektóre w ciągłych luksusach (gdyż nie są to rozkłady normalne). W szczególności widać to w GDP, gdzie większość krajów znajduje się bardzo nisko, a tylko elitarnych kilka osiąga duże wartości.

7. Korelacje między zmiennymi

Do porównania korelacji atrybutów został zastosowany współczynnik korelacji Pearsona. Na wykresie poniżej przedstawione są (razem z odbiciem lustrzanym) wszystkie wartości. Warto wspomnieć, iż krzyżyki zaznaczone są tam, gdzie korelacja jest mniejsza niż 0.05 (tz. poziom istotności).

Możemy zauważyć, iż liczba populacji nie wpływa na długość życia. Także rok oraz procent rządowych wydatków nie jest bardzo mało skorelowany z długością życia. Może to być wynikiem różnych krajów “rozwiniętych”, gdzie procent ten jest utrzymywany na stopniu, który zapewnia odpowiedni poziom opieki medycznej (procent może się różnić, ale opieka pozostaje podobna).

library(ggcorrplot)
library(ggplot2)

corr <- round(cor(numericData, method = "pearson"), 1)
p.mat <- cor_pmat(numericData)

dev.new(width=8, height=8)

gg <- ggcorrplot(corr, p.mat = p.mat, hc.order=TRUE,
                 type='full', tl.cex=8,
                 sig.level = 0.05, insig= "blank")

fig <- ggplotly(gg)
fig

8. Długość życia a kraj

Średnia długość życia dla różnych krajów została przedstawiona na wykresie słupkowym poniżej. Kraj o statusie “Developing” zostały oznaczone na niebiesko, podczas gdy “Developed” – na zielono. Możemy zauważyć, iż kraje rozwinięte są przeważnie w czołówce średniej długości życia, nie znaczy to jednak, iż kraje rozwijające się są gorsze – niektóre z nich, takie jak Kanada i Francja są lepsze od większości krajów rozwiniętych.

library(plotly)

countryNames <- unique(cleanedData$Country)
countryData <- cleanedData %>%
  filter(Life.Expectancy != 0) %>%
  group_by(Country, Status) %>%
  dplyr::summarize(Mean.Life.Expectancy = mean(Life.Expectancy), .groups = 'drop')

cbPalette <- countryData$Status
cbPalette[cbPalette == "Developed"] <- "green"
cbPalette[cbPalette == "Developing"] <- "#56B4E9"

fig <- plot_ly(
  countryData,
  x = ~Country,
  y = ~Mean.Life.Expectancy,
  type = 'bar',
  marker = list(color = cbPalette)
  )

fig

9. Regresor

Wytrenowałem regresor – Random Forest Regressor – na wszystkich danych numerycznych tego zbioru. Została zastosowana 2-krotna walidacja krzyżowa z 5-cio krotnym powtórzeniem. Dane zostały przeskalowane.

library(caret)

setConsistency()
inTraining <-
    createDataPartition(
        y = numericData$Life.Expectancy,
        p = .75,
        list = FALSE,)

trainingSet <- numericData[inTraining,]
testingSet <- numericData[-inTraining,]

rfGrid <- expand.grid(mtry = seq(2, ncol(numericData) - 1))
gridCtrl <- trainControl(
    method = "repeatedcv",
    number = 2,
    repeats = 5)

fitTune <- train(Life.Expectancy ~ .,
             data = trainingSet,
             method = "rf",
             metric = "RMSE",
             preProc = c("center", "scale"),
             trControl = gridCtrl,
             tuneGrid = rfGrid,
             ntree = 10)
fitTune
## Random Forest 
## 
## 1334 samples
##   19 predictor
## 
## Pre-processing: centered (19), scaled (19) 
## Resampling: Cross-Validated (2 fold, repeated 5 times) 
## Summary of sample sizes: 666, 668, 667, 667, 667, 667, ... 
## Resampling results across tuning parameters:
## 
##   mtry  RMSE      Rsquared   MAE     
##    2    2.865870  0.8981423  2.008806
##    3    2.611478  0.9138562  1.823620
##    4    2.496997  0.9209252  1.729448
##    5    2.414434  0.9260085  1.651294
##    6    2.334129  0.9300500  1.604995
##    7    2.349975  0.9291291  1.599798
##    8    2.281503  0.9329390  1.568389
##    9    2.303874  0.9315994  1.572533
##   10    2.300910  0.9315208  1.563833
##   11    2.288248  0.9320849  1.549882
##   12    2.274257  0.9330773  1.553570
##   13    2.288128  0.9322114  1.554789
##   14    2.282676  0.9322473  1.545455
##   15    2.297981  0.9313793  1.561153
##   16    2.286173  0.9322201  1.551356
##   17    2.328094  0.9294172  1.578184
##   18    2.300977  0.9314149  1.557778
##   19    2.324170  0.9297220  1.572661
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 12.
library(ggplot2)
library(plotly)
gg <- ggplot(fitTune) + theme_bw()
fig <- ggplotly(gg)
fig

Jak widać, regresor zwraca bardzo dobre wyniki, na poziomie 0.95 R2. Oznacza to prawdopodobnie, iż średnia jakość życia może zostać bardzo prosto oszacowana na podstawie danych o kraju.

rfTuneClasses <- predict(fitTune,
                         newdata = testingSet)

cat("\nTesting set:\n")
## 
## Testing set:
postResample(
  as.numeric(rfTuneClasses),
  testingSet$Life.Expectancy)
##      RMSE  Rsquared       MAE 
## 1.9878760 0.9487156 1.2810903

10. Najważniejsze atrybuty

Jest to jednakże klasyfikator, które prawdobnie jest bardzo szczegółowy, gdyż ilość wybranych predyktorów wynosi 12. Kolejny regresor jest dokładnie taki sam jak poprzedni, aczkolwiek dane do jego treningu do 6 najbardziej korelujących z długością życia atrybutów.

library(caret)
library(dplyr)

bestCorr <- cor(numericData) %>%
  as.data.frame() %>%
  dplyr::select(Life.Expectancy) %>%
  abs(.) %>%
  arrange(desc(Life.Expectancy)) %>%
  slice(1:7)

lessData <- numericData %>%
  dplyr::select(as.vector(row.names(bestCorr)))

setConsistency()
inTraining <-
    createDataPartition(
        y = lessData$Life.Expectancy,
        p = .75,
        list = FALSE,)

trainingSet <- lessData[inTraining,]
testingSet <- lessData[-inTraining,]

rfGrid <- expand.grid(mtry = seq(2, ncol(lessData) - 1))
gridCtrl <- trainControl(
    method = "repeatedcv",
    number = 2,
    repeats = 5)

fitLessTune <- train(Life.Expectancy ~ .,
             data = trainingSet,
             method = "rf",
             metric = "RMSE",
             preProc = c("center", "scale"),
             trControl = gridCtrl,
             tuneGrid = rfGrid,
             ntree = 10)
fitLessTune
## Random Forest 
## 
## 1334 samples
##    6 predictor
## 
## Pre-processing: centered (6), scaled (6) 
## Resampling: Cross-Validated (2 fold, repeated 5 times) 
## Summary of sample sizes: 666, 668, 667, 667, 667, 667, ... 
## Resampling results across tuning parameters:
## 
##   mtry  RMSE      Rsquared   MAE     
##   2     2.304066  0.9315043  1.587042
##   3     2.263488  0.9334297  1.547365
##   4     2.254369  0.9340318  1.517265
##   5     2.278465  0.9325793  1.527133
##   6     2.291570  0.9318532  1.543398
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 4.
rfLessTuneClasses <- predict(fitLessTune,
                         newdata = testingSet)

cat("\nTesting set:\n")
## 
## Testing set:
postResample(
  as.numeric(rfLessTuneClasses),
  testingSet$Life.Expectancy)
##      RMSE  Rsquared       MAE 
## 1.9010933 0.9529699 1.1967020

Jak możemy zauważyć, liczba predyktorów zmniejszyła się do 4, a wynik wzrósł. Nie jest to wzrost statystycznie coś mówiący, jednakże sam fakt, iż wynik pozostał na tym samym poziomie 0.95 R2 mówi nam, iż ten klasyfikator powinien być dużo lepszy w generalizacji dla nowych danych.

library(ggplot2)
library(plotly)
gg <- ggplot(fitLessTune) + theme_bw()
fig <- ggplotly(gg)
fig